%% process_mesh
function mesh = process_mesh(pde, varargin)
if ~isempty(varargin)
    h = varargin{1};
    if h > (pde.domain(2)-pde.domain(1)) || h > (pde.domain(4)-pde.domain(3))
        error("h is too large.");
    end
    mesh.Nx = (pde.domain(2)-pde.domain(1))/h;
    mesh.Ny = (pde.domain(4)-pde.domain(3))/h;
    [mesh.node, mesh.elem] = process_square(pde.domain, mesh.Nx, mesh.Ny);
    [mesh.edge, mesh.bde, mesh.elem_ed] = process_edge(mesh);
    mesh.Np = size(mesh.node,1);
    mesh.Ne = size(mesh.elem,1);
    mesh.Ned = size(mesh.edge,1);
    mesh.Nbd = 4;
else
    fileID = fopen(pde.meshfile, "r");
    data = textscan(fileID, "%f %f %f %f", "Delimiter", " ", "MultipleDelimsAsOne", true);
    fclose(fileID);
    Np = data{1}(1);
    Ne = data{2}(1);
    mesh.node = [data{1}(2:Np+1), data{2}(2:Np+1)];
    mesh.elem = [data{1}(Np+2:Np+Ne+1), data{2}(Np+2:Np+Ne+1), data{3}(Np+2:Np+Ne+1)];
    mesh.bdary = [data{1}(Np+Ne+2:end), data{2}(Np+Ne+2:end), data{3}(Np+Ne+2:end)];
    [mesh.edge, mesh.bde, mesh.elem_ed] = process_edge(mesh);
    mesh.Np = Np;
    mesh.Ne = Ne;
    mesh.Ned = size(mesh.edge,1);
    mesh.Nbd = max(unique(mesh.bdary(:,3)));
end
end

%% process_square
function [node, elem] = process_square(domain, Nx, Ny)
[X, Y] = meshgrid(linspace(domain(1), domain(2), Nx+1), linspace(domain(3), domain(4), Ny+1));
node = [X(:), Y(:)];
Np = size(node,1);
k = (1:Np-Ny-1)';
cut = (Ny+1)*(1:Nx);
k(cut) = [];
elem = [k, k+Ny+1, k+1;
    k+1, k+Ny+1, k+Ny+2];
elem = sortrows(elem);
end

%% process_edge
function [edge, new_bde, elem_ed] = process_edge(mesh)
elem = mesh.elem;
Ne = size(elem,1);
% Np = size(mesh.node,1);
% Ned = Ne+Np-1;
% edge = zeros(Ned,2);
% elem_ed = zeros(Ne,3);
% edge_el = edge;
% edgeMap = containers.Map("KeyType", "char", "ValueType", "any");
% count = 0;
% t = [1 2;
%     2 3;
%     3 1];
% for e = 1:Ne
%     tempe = reshape(elem(e,t), 3, 2);
%     for m = 1:3
%         key = sprintf("%d,%d", tempe(m,2), tempe(m,1));
%         if isKey(edgeMap, key)
%             idx = edgeMap(key);
%             elem_ed(e,m) = -idx;
%             edge_el(idx,2) = e;
%         else
%             count = count + 1;
%             edge(count,:) = tempe(m,:);
%             elem_ed(e,m) = count;
%             edge_el(count,:) = [e,0];
%             new_key = sprintf("%d,%d", tempe(m,1), tempe(m,2));
%             edgeMap(new_key) = count;
%         end
%     end
% end
% bde = [edge_el(edge_el(:,2)==0,1), edge(edge_el(:,2)==0,:), find(edge_el(:,2)==0)];
edge_all = zeros(3*Ne,3);
t = [1 2;
    2 3;
    3 1];
k = 1;
for e = 1:Ne
    edge_all(k:k+2,:) = [e*ones(3,1), reshape(elem(e,t), 3, 2)];
    k = k+3;
end
[~, ib, ic] = unique(sort(edge_all(:,2:3),2), "rows", "stable");
k = accumarray(ic,1);
idx = ic;
idx((k(ic)==2)) = -idx((k(ic)==2));
idx(ib) = ic(ib);
elem_ed = reshape(idx, 3, [])';
edge = edge_all(ib,2:3);
bde = [edge_all(ib(k==1),:), find(k==1)];
bdeMap = containers.Map("KeyType", "double", "ValueType", "any");
for i = 1:size(bde,1)
    bdeMap(bde(i,2)) = i;
end
if isfield(mesh, "bdary")
    bdary = mesh.bdary;
    new_bde = zeros(size(bde));
    for i = 1:size(bde,1)
        new_bde(i,:) = bde(bdeMap(bdary(i,1)),:);
    end
    new_bde = [new_bde, bdary(:,3)];
else
    new_bde = zeros(size(bde));
    new_bde(1,:) = bde(1,:);
    for i = 2:size(new_bde,1)
        idx = bdeMap(new_bde(i-1,3));
        new_bde(i,:) = bde(idx,:);
    end
end
end